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The recently developed tensor renormalization-group (TRG) method provides a highly precise 
technique for deriving thermodynamic and critical properties of lattice Hamiltonians. The TRG is 
a local coarse-graining transformation, with the elements of the tensor at each lattice site playing 
the part of the interactions that undergo the renormalization-group flows. These tensor flows are 
directly related to the phase diagram structure of the infinite system, with each phase flowing to a 
distinct surface of fixed points. Fixed-point analysis and summation along the flows give the critical 
exponents, as well as thermodynamic functions along the entire temperature range. Thus, for the 
ferromagnetic triangular lattice Ising model, the free energy is calculated to better than 10 -5 along 
the entire temperature range. Unlike previous position-space renormalization-group methods, the 
truncation (of the tensor index range D) in this general method converges under straightforward 
and systematic improvements. Our best results are easily obtained with D = 24, corresponding to 
4624-dimensional renormalization-group flows. 
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I. INTRODUCTION 

The tensor renormalization-group (TRG) method, re- 
cently introduced by Levin and Nave [l[ is versatile, ac- 
curate, and conceptually interesting. The versatility of 
TRG comes from being a genuinely local renormalization- 
group transformation, that is a mapping between lo- 
cal Hamiltonians on the original and coarse-grained lat- 
tices, expressed in terms of tensors at sites of each lat- 
tice. Thus, although TRG has been demonstrated by 
obtaining very accurate phase transition temperatures 
and magnetization curves for frustrated and unfrustrated 
classical Ising models on the triangular lattice, we be- 
lieve that it can be extended to yield the critical prop- 
erties and entire thermodynamics of complicated multi- 
critical systems, including quenched random systems, by 
renormalization-group flow analysis. 

In the present work, we show that the renormalization- 
group flows of the tensor elements can be analyzed in 
the same manner as the flows of Hamiltonian interaction 
parameters, which was not done in the previous work. 
Each thermodynamic phase region corresponds to a basin 
of attraction in the space of tensor elements, and criti- 
cal temperatures and exponents (not obtained before) 
are deduced from characteristics of the boundaries be- 
tween the different basins. Practically, this means that 
the phase diagram, thermodynamic, and critical proper- 
ties of the infinite system can be derived directly from 
the flows. Moreover, the flows themselves have an in- 
teresting and unconventional nature, since a phase re- 
gion does not flow to single isolated sink, but rather to 
a continuous surface of fixed points. Using the ferro- 
magnetic triangular lattice Ising model as an example, 
we show that the accuracy of the calculated free energy, 
critical temperature, and thermal critical exponent can 
be easily and systematically improved on by increasing 



D, the cutoff on the index range of the tensors. In this 
ability to converge toward exact values with larger cut- 
offs, TRG is more general and dependable than position- 
space renormalization-group techniques @, H, 0] whose 
successes have been based on system-specific heuristics 
[E B • The TRG approach combines the straightfor- 
ward interpretative framework of traditional renormal- 
ization group analysis of flows in a parameter space — 
with the accuracy of techniques depending on finite-size 
scaling of large systems. 



II. TENSOR RENORMALIZATION-GROUP 
TRANSFORMATION 

The TRG method [l] can be applied to any classical 
lattice Hamiltonian satisfying the following conditions: 
(1) it can be expressed in terms of degrees of freedom on 
the bonds of the lattice; (2) the Boltzmann weight of a 
configuration can be written as the product of individual 
Boltzmann weights for each lattice site, that only depend 
on the bond variables adjoining the site. This is a broad 
category including many common statistical physics sys- 
tems like the Ising and Potts models, both of which can 
be mapped onto this form through duality transforma- 
tions, as well as all vertex models. For such Hamiltoni- 
ans, the partition function can be written as a tensor net- 
work [1, 0] : We start with a lattice of N sites where each 
site has coordination number q and each bond is in one of 
d possible states. The Boltzmann weight of an individual 
site depends on the configuration of the q bonds meeting 
at the site and can be written as a tensor Tj 1 j 2 ...j , where 
each index i a runs from 1 to d. For configurations of bond 
variables that are not allowed, the corresponding clement 
of T is zero. The tensor is real- valued and cyclically 
symmetric. The transformation described below works 
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FIG. 1: (a) Graphical representation of the tensors, (b) 
Rewiring, (c) Decimation, (d) Renormalization group trans- 
formation applied to the entire hexagonal lattice, with the 
first arrow showing the rewiring step, and the second arrow 
showing the decimation. 



also for the more general case of complex-valued, cycli- 
cally symmetric tensors, onto which our original tensor 
network system is mapped after a single renormalization 
step. The partition function is a product over the N site 
tensors, with each index contracted between two different 
tensors (since each bond is shared between two sites), 

d 

Z = ^ t Ti 1 i 2 ...i q Ti 1 i r ...i B Ti 2 i t ...i u ■ ■ ■ . (1) 

•lv,!M = l 

The position-space renormalization-group transforma- 
tion of this system allows us to express the partition 
function equivalently as a tensor network over a coarse- 
grained lattice with N' = N/b 2 vertices, in two spatial 
dimensions, with a length rescaling factor b > 1. It is 
accomplished in two steps, which we call rewiring and 
decimation. For simplicity, we shall focus in this work 
on applying the transformation to the hexagonal lattice, 
though this two-step procedure can easily be adapted 
to a variety of two-dimensional lattices, including the 
square [l[ and Kagome lattices. 

Rewiring: Graphically, let us represent the tensor Tijk 
at each point of the hexagonal lattice as a three-legged 
vertex as shown on the left of Fig. [2(a), and a contraction 
of an index between two tensors as a connection between 
two vertex legs. We use the convention that the order 
of indices in the tensor matches the counterclockwise or- 
dering of the labels on the vertex legs. The rewiring step 
for a pair of neighboring tensors consists of reconnecting 
the bonds in the manner of Fig. (Hb), rewriting them as 
a contraction of two new tensors S, 

d d 2 

^ ^ T-ijk^-kUn ^ ^ SyniuSjiu . (2) 
fc=l v=l 



The first two indices in the tensor S m i V run up to d, but 
the third index v runs up to d 2 . We shall distinguish 
these types of indices by a Greek letter label, and graph- 
ically depict them as thick legs [Fig. [2(a) center]. To see 
that such a rewiring is possible, let us introduce compos- 
ite indices a = (m, i) and (3 = (j, I), and write the tensor 
contraction on the left-hand side of Eq. @ as a d 2 x d 2 
matrix M a p = J^kTijkTkim- Then Eq. @ becomes 

Ma/3 = 2^ SctvSpu) (3) 
v 

or M — SS T . We can find S using the fact that 
M = M T , as can be checked using the cyclical sym- 
metry of the T tensors. Any symmetric matrix M ad- 
mits a variant of singular value decomposition known as 
Takagi factorization [l(|: M = UYJJ 1 ', where U is a 
d 2 x d 2 unitary square matrix and £ is a d 2 x d 2 diag- 
onal matrix containing the singular values of M . If £„„ 
is the vth singular value (assumed ordered from largest 
to smallest with increasing v), then the elements of S 
are given by S av = \fT^JJ av . The resulting decomposi- 
tion M = SS T is not unique, because we can replace S by 
SO, where O is a complex orthogonal matrix (OO t = I). 
However, we shall always use the S given directly by the 
the Takagi factorization, so if the vth singular value Yi uv 
is nondegenerate, the z^th column of S is uniquely deter- 
mined up to a factor of ±1. The rewiring procedure de- 
scribed above is applied globally to the entire hexagonal 
lattice as shown by the first arrow in Fig. HKd), grouping 
the T tensors into pairs, and replacing each pair by S 
tensors. 

Decimation: The second step in the renormalization 
procedure consists of tracing over the degrees of free- 
dom in the triangular clusters that are formed after the 
rewiring. As illustrated in Fig. [TJc), each such cluster 
can be replaced by a single renormalized tensor T', 

d 

^ Shm^S m iySih8 = T^fi u . (4) 
m,i,/i— 1 

Introducing the notation for the d x d matrix with 
elements Shmj at fixed 7, then the above equation has 
the form 

Tr[5 (7) 5 (l/) S* (,5) ] =T' lSv . (5) 

It is clear from Eq. ([5]) that the renormalized tensor will 
be cyclically symmetric. Since the 5* tensor is constructed 
from the unitary matrix U, the elements of T' will in 
general be complex. The second arrow of Fig. [2(d) shows 
the decimation applied to every triangular cluster, and 
the result is a hexagonal lattice of T" tensors with the 
lattice spacing larger by a factor of b = -\/3 compared to 
the original hexagonal lattice. 

The renormalization-group transformation described 
so far is an exact mapping of one tensor network onto an- 
other, preserving the partition function Z . However, the 
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renormalized tensors have a more complicated structure 
than the original ones, since the indices of T' &v run up to 
d 2 . If this procedure is repeated, the index range grows 
exponentially with iteration number, making numerical 
implementation impractical. This is analogous to the dif- 
ficulty encountered when applying naive position-space 
renormalization to Hamiltonians on lattices, where the 
number of couplings in the Hamiltonian grows with each 
iteration. Some kind of approximate truncation is re- 
quired to keep the complexity of the renormalized model 
bounded. In our case this truncation can be done in 
a straightforward and systematic fashion by setting an 
upper bound D for the index range. Rather than use 
the full d 2 x d 2 matrix S al/ to calculate T", we use only 
the first d columns, where d = min (d 2 ,D). Using this 
truncated d 2 x d matrix, which we call S, means that 
the rewiring step is implemented only approximately, 
M « SS T . However, this is the best approximation pos- 
sible, since the first d columns of S correspond to the d 
largest singular values of M. The decimation step of Eq. 
(d|) is still carried out exactly, but with S replaced by S. 
The resulting tensor T' Sll has indices that run up to d. 

III. RENORMALIZATION-GROUP FLOWS 
AND THERMODYNAMIC BEHAVIOR OF THE 
TRIANGULAR LATTICE ISING MODEL 

In order to illustrate the nature of the renormalization- 
group flows resulting from the transformation outlined 
above and the methods by which thermodynamic infor- 
mation can be extracted from them, we turn now to a spe- 
cific example: the triangular lattice ferromagnetic Ising 
model. This model is mapped onto a tensor network 
through a duality transformation, giving an equivalent 
partition function in terms of bond variables on a hexag- 
onal lattice, [lj Each bond in this dual lattice has two 
states, ctj. = 1 or era = — 1, corresponding to a bond be- 
tween parallel or antiparallel Ising spins. These states 
respectively contribute energies — J < and J to the to- 
tal Hamiltonian 7i. The tensor for each lattice site has 
the form 

^ = ^^'-^(^ + 1), (6) 

where j3 = l/fc^T, and the factor multiplying the ex- 
ponential is a projection operator that is equal to 1 for 
allowed configurations of the bonds and otherwise. For 
simplicity we set J/ks = 1, effectively measuring tem- 
peratures in units of J /kg. If the triangular lattice has 
N sites, the dual hexagonal lattice has 2N sites, and the 
partition function is a contraction of the 2N site tensors. 

Iterating the renormalization-group transformation, 
we consider renormalization-group flows in the space of 
tensor element amplitudes |T a/ 3 7 |. For a cyclically sym- 
metric tensor Ty-j. with index range d, the maximum 
number of distinct elements is d(2 + d 2 )/3. After the 
first few iterations, the index range of the renormalized 



tensors reaches the cutoff D and remains there, so that 
the flows in the subsequent steps are in a space with di- 
mension D(2 + D 2 )/3. In this work we look at cutoffs in 
the range D — 4 to 24, so the most complex flow space 
that we investigate is 4624-dimensional. A large number 
of the D(2 + D 2 )/3 tensor elements are zero, correspond- 
ing to disallowed bond configurations, typically making 
up between 50-55% of the total in the later steps. The 
locations of the zero elements usually stay the same for 
a few steps at a time, but can shift abruptly during the 
flow. To minimize roundoff errors which can affect the 
results at later iterations, the singular value decomposi- 
tion and matrix manipulations involved in TRG are all 
implemented at quadruple precision, using a customized 
version of the LAPACK library [HI- 

Due to the contraction during the decimation step, the 
tensor elements tend to increase in magnitude rapidly 
with each iteration, which can lead to numerical prob- 
lems. These are avoided similarly to the procedure 
used in global strong-coupling flows in position-space 
renormalization-group theory. The physics described by 
the tensor network is unaffected by an overall constant 
multiplying each element in the tensor. Thus at each step 
we factor out the absolute value of a nonzero element, 
for example T — \Tm\T, yielding a reduced tensor T. 
The Tm element is a convenient choice, since it remains 
nonzero throughout the renormalization-group flow. The 
renormalization-group transformation is then applied to 
the reduced tensor T, giving a renormalized tensor T'. 
Let us denote T as and T as TW. The factorization 
and renormalization-group transformation are iterated, 
so at the nth step we have a tensor T^> = \T^l |T» . As 
we shall see below, with this added factorization step the 
elements of T^ n > tend to finite limiting values as n — ► oo. 

Additionally, keeping track of the factors \T^l \ during 
the flow, as in standard position-space renormalization- 
group theory, allows us to easily determine the free en- 
ergy of the system in the thermodynamic limit. [l2| Defin- 
ing G^™) = In | T^"^ | , the partition function is a con- 
traction of 2N tensors T^°\ which we write schemat- 
ically as Z = (T^) 2N = e 2NGW {f(°)) 2N . After a 
single renormalization-group step, this becomes Z = 

e 2ATG<°>( T (l))2A73 = e 2JV(G»+GW/3)(f(l))2W/3 i Thus af _ 

ter n steps we find 

Z = e 2A^r=o3- < G« (f(n)y N/3 " (?) 

The free energy per original site j3f = —N^ 1 In Z is then 
^/ = -2^3- i GW--ln(fW) . (8) 

i=0 

As n — > oo the elements of go to finite fixed val- 
ues, and thus for large N the second term on the right 
in Eq. ((5|) becomes negligible. So the final result for 
the free energy is 0f — — 2^^ 3 _l GW. In practice 
the transformation is iterated for n steps until the value 
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FIG. 2: The relative error of / (D) (T), the free energy per site 
of the triangular lattice Ising model calculated using TRG 
at various cutoffs D = 4 to 24, compared to the exact free 
energy /exact (T), as a function of temperature T. The dotted 
line shows the exact critical temperature T c = 4/ In 3. 



of GW has converged to the limit G^ 00 ) within the de- 
sired numerical precision. Then (3f sa — 2 Y^i=o 3~*G^ — 
2G (n) Y,rL„ = ~ 2 Y%=o 3 ~ l G (4) ~ 3 1_n G (n) . 

Even for small cutoffs D, the results of this free energy 
calculation can be remarkably accurate, as seen in Fig. 
[2j which shows the relative error of f (D) (T), the free 
energy per site using cutoff D, compared to the exact 
value /exact (T) . The latter is calculated by numerical 
evaluation of the integral solution in Ref. [13[. Already 
at D = 4 the TRG free energy is within 0.09% of the 
exact value at all temperatures, and is considerably more 
accurate than this away from the critical region near T c — 
4/ In 3. As noted in Ref. [l[ , the TRG method is expected 
to behave worst at criticality, and this is indeed what 
we see for all D, with the relative error curves peaked 
near T c . Yet even here there is significant improvement 
as we go to larger cutoffs. At D = 24, the largest cutoff 
examined, the maximum error is 0.0007%. Overall, going 
from D = 4 to D = 24 we get an improvement between 
two and three orders of magnitude in the precision of the 
free energy result. 

Further thermodynamic information can be gleaned by 
looking in detail at the behavior of the renormalization- 
group flows. The elements \T$\ go to finite fixed 
values \T*j k \ as n — > oo, but unlike typical position- 
space renormalization-group transformations, there are 
no unique fixed points which act as a basins of attraction 
for the low- and high-temperature phases. Instead, at 
each different temperature T the system flows to a dif- 
ferent fixed point \T*- k (T)\. However there is a way to 
distinguish the low- and high-temperature phases, since 
each flows to a different continuous surface of fixed points 
in the space of tensor elements. There is a boundary tem- 
perature Tc°\ depending on the cutoff D, such that for 
T > T C (D) the flows tend to one fixed surface, while for 
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FIG. 3: The relative error of the critical temperature T c 
for the triangular lattice Ising model calculated from TRG 
flows at various cutoffs D — 4 to 24, compared to the exact 
T c = 4/ In 3. 



T < Tc they tend to the other. As D becomes larger, 
Tc D ^ gives a rapidly converging estimate of the exact crit- 
ical temperature T c = 4/ In 3. Fig. [3] shows the relative 
error of compared to T c for various D. While the 
decrease in error is not monotonic with D, there is an 
overall trend which appears to be roughly exponential in 
the cutoff, so that the error at D = 24, 0.02%, is two 
orders of magnitude smaller than at D = 4. 

This segregation of the flows between two distinct fixed 
surfaces can be seen directly in Fig. [4] which plots one 
tensor element |T^| as a function of iteration number n 
for D = 8 and D — 12. The flows are shown for different 
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are two curves for each |AT|, one corresponding to T > 
T C (D) , and the other to T < T C (D) . The two curves stay 
close to one another for a number of iterations, but then 
veer off in opposite directions and reach different fixed 
values |Xx*i|. These fixed values change continuously as 
| AT | is varied, as they map out a slice of the two fixed 
surfaces corresponding to the low- and high-temperature 
phases. 

As |AT| gets smaller, there is an interesting difference 
in the flow behaviors of the D = 8 and D = 12 cases. 
For D = 8, as seen in Fig. [DJa) and in more detail in 
Fig. 0Jb), the AT > and AT < curves are nearly 
horizontal before turning away to their respective fixed 
surfaces. This is true in fact for all the |T Q( g 7 | elements, 
since the flows are attracted to a unique critical fixed 
point, which we denote |T^ |. The smaller the value 
of | AT |, the larger the number of iterations which are 
spent in the vicinity of |T°oJ before flowing to one of 
the fixed surfaces. Once we are in the vicinity of the crit- 
ical fixed point, we can isolate it to high precision using 
a Newton-Raphson procedure. The analysis of the fixed 
point proceeds just as in a standard renormalization- 
group approach: we calculate a recursion matrix, whose 
eigenvalues can be related to the critical exponents of the 
system. To do this, let us denote the nonzero elements of 
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FIG. 4: The behavior of tensor element |T ln | as a function 
of iteration number n. The different data sets correspond to 
flows at different temperatures T near , with each symbol 
in the legend corresponding to a value of [AT| = \T — Tc D ^\ 
between 10 -2 and 10 -10 . For each |AT| there are two data 
sets, one for AT > connected by solid lines, and the other for 
AT < connected by dashed lines. For clarity, the flows are 
only shown up to the iteration where they have approximately 
converged to a fixed value, (a) IT^-q | flows with cutoff D = 8. 
(b) Same as in (a), but zoomed in to see more clearly the 
flows for AT < 0. (c) [T^l flows with cutoff D = 12. 



\T£pj\, not related by cyclical symmetry, as K\ through 
K m . In the case of D = 8, m = 80. For small pertur- 
bations away from the critical fixed point, the number 
and locations of these nonzero elements stay the same 
after a renormalization-group transformation, which al- 
lows us to numerically evaluate the mx m recursion ma- 
trix Rij = dK[jdKj. Writing the eigenvalues of R in 
the form b Vi , i — 1, . . . , m, we find only one eigenvalue 
where yi > 0, as expected at a critical fixed point. This 
relevant eigenvalue, which we denote j/t, is related to the 
specific heat critical exponent a by a = (2yr — 2)/yr- 
For D = 8, y T = 1.01543, giving a = 0.03039, which 
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FIG. 5: Specific heat per site C of the triangular lattice Ising 
model, calculated using TRG for various cutoffs D — 4 to 
24 as a function of T — Ti D ' in the high-temperature phase 
near criticality. The gray curves superimposed on the data 
points are best-fit curves of the form A + B(T — Ti 15 ') -0 , 
with parameters A, B, and a. The thick black line is the 
exact specific heat as a function of T — T c . 



D 


T (n) 


T C (D) -T C 


yr 


4 


3.73840 


9.7 x 10 -2 


1.01543 


8 


3.60873 


-3.2 x 10~ 2 


1.01543 


12 


3.64958 


8.6 x 10 -3 


1.0644 ± 0.0048 


16 


3.63847 


-2.5 x 10~ 3 


1.0574 ± 0.0022 


20 


3.64147 


5.1 x 10 -4 


1.0165 ±0.0032 


24 


3.64183 


8.7 x 10" 4 


1.0089 ±0.0075 


exact 


3.64096 
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TABLE I: Critical properties calculated using the TRG ap- 
proach, at various cutoffs D, compared to the exact values in 
the last row. The thermal eigenvalues yr for D — 4 and 8 are 
calculated from the recursion matrix evaluated at the critical 
fixed point. For D > 8 the estimate for j/t is from the best-fit 



result to the specific heat near Ti , as plotted in Fig. [S] 



compare well to the exact values of j/t = 1 and a = 0. 
We can check this result using an alternative approach, 
by calculating the specific heat from derivatives of the 
calculated free energy for small T - T^ D \ and we find 
that the singularity in the specific heat agrees with the a 
derived from the thermal eigenvalue yx of the recursion 
matrix. 

For the six values of the cutoff D where we investigated 
the near-critical flows in detail, D = 4, 8, 12, 16, 20, 24, 
we were able to isolate the critical fixed point for D = 4 
and 8, namely in 24- and 196-dimensional flow spaces. 
Despite the different dimensionalities of the flow spaces, 
both cases yielded the same eigenvalue yT, within the 
precision of 5 decimal places. Higher values of D showed 
very different flow behaviors, as exemplified in Fig. HJc) 
for D = 12. Here the AT > and AT < curves 
for |T^"j| do not stay nearly horizontal before diverging 
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to the fixed surfaces: regardless of how small |AT| is 
made, the flows do not gravitate toward a unique critical 
fixed point, but map out a continuous spectrum of points 
which attract the flows before the two curves spread out. 
Using an arbitrary precision version of the TRG algo- 
rithm implemented in Mathematical we checked |AT| 
values as small as 10~ 43 without finding convergence to- 
ward a critical fixed point. Nevertheless, in these cases 
for D > 8 we can still extract the critical behavior, by 
resorting to the alternative approach mentioned earlier: 
looking directly at the specific heat C per site near Tj D ^ 
as obtained from the numerical derivative of the calcu- 
lated free energy. We plot C within the high-temperature 
phase for various D in Fig. El for T - T^ D) between 1CT 7 
and 10 -4 . For each D the data points are fit to a function 
A + B(T — Tc D ^)^ a , which provides an accurate descrip- 
tion of the singularity, and the best-fit value of a is used 
to determine the thermal eigenvalue y-p. The results are 
listed in Table [J For D = 12, y T = 1.0644 ± 0.0048, 
worse than our fixed-point determined values at D = 4 
or 8, but this value improves as D is increased, reaching 
y T = 1.0089 ± 0.0075 at D = 24. Moreover, as seen in 
Fig. O our calculated C curve for D = 24 nearly overlaps 
the exact C as a function of T — T Cl which diverges with 
a logarithmic singularity at the critical point. 



the triangular lattice Ising model, the low- and high- 
temperature phase regions are basins of attraction for 
two distinct surfaces of fixed points. The boundary be- 
tween these basins defines a critical temperature Tc D \ 
dependent on the TRG cutoff D. At small cutoffs such 
as D = 4 and 8, the flows near the boundary between 
the basins are controlled by a critical fixed point, while 
at higher D the flows show more complicated behavior, 
never converging at a unique point. In the former case 
the thermal exponent yx is found from the eigenvalues of 
the recursion matrix at the critical fixed point, while in 
the latter we can deduce the exponent from the scaling of 
the calculated specific heat near T^ D) . The free energy 
at all temperatures systematically converges to the ex- 
act Ising result with increasing D, particularly fast away 
from the critical region. For the critical properties the 
improvement is not monotonic in D, but both the criti- 
cal temperature and the exponent yx tend toward 
the exact Ising values at larger cutoffs. With very modest 
computational effort, the TRG method provides an accu- 
rate portrait of global phase diagram characteristics. It 
thus warrants further study, both on applications to other 
two-dimensional lattice models, and possible generaliza- 
tion to systems with quenched randomness [14| and/or 
higher spatial dimensions. 



IV. CONCLUSION 
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